﻿#include <iostream>
#include <iomanip>
#include <vector>
#include <Eigen/Sparse>
#include <Eigen/SparseLU>
#include <fstream>


using namespace std;

enum BASIS_TYPE//basis function type
{
	ONE_DIM_LINEAR,//1D linear
	ONE_DIM_QUADRATIC//1D quadratic
};

enum BOUNDARY_TYPE//boundary type
{
	Dirichlet,
	Neumann,
	Robin,
	Other
};

struct boundary_condition
{
	BOUNDARY_TYPE type;//define the boundary type
	int glo_node_index;//define the global index of the boundary node
	vector<double> value;//define the boundary value, it can be very complex and define by a function
	vector<double> outer_normal_vector;//outer normal vector of the boundary used by Robin and Neumann conditions
};

typedef Eigen::SparseMatrix<double> SpMat;//declares a sparse matrix type of double
typedef Eigen::VectorXd Rsv;//declares a right side value type
//typedef Eigen::Triplet<double> T;//declares a triplet to stores the matrix index and values, then convert it to sparse matrix(with high performance)


//mesh generate and get the T, P matrix
//P and T for mesh nodes
//Pb_test and Tb_test for finite elements nodes corresponding to basis test functions
//Pb_trial and Tb_trial for finite elements nodes corresponding to basis trial functions
//For the nodal finite element basis functions, P=Pb, T=Tb
//In the Gherkin method, Pb_test=Pb_trial and Tb_test=Tb_trial
//For general purpose, handle the three type of matrix respectively
//return: P,T,Pb_trial,Tb_trial,number_of_local_basis_trial,Pb_test,Tb_test,number_of_local_basis_test
bool mesh_generate(BASIS_TYPE basis_type_trial, BASIS_TYPE basis_type_test,double x_min, double x_max, int number_of_elements,
	vector<vector<double>> &P, vector<vector<int>> &T,
	vector<vector<double>> &Pb_trial,vector<vector<int>> &Tb_trial, int &number_of_local_basis_trial,
	vector<vector<double>> &Pb_test, vector<vector<int>> &Tb_test, int &number_of_local_basis_test)
{
	if (x_min >= x_max)
	{
		cout << "Input x_min and x_man error, please check" << endl;
		return false;
	}

	//P and T matrix
	P.push_back(vector<double>{x_min});
	for (int i = 0; i < number_of_elements; ++i)
	{

		P.push_back(vector<double>{(x_max - x_min) / number_of_elements * (i + 1)});
		vector<int> temp;
		temp.push_back(i);
		temp.push_back(i + 1);
		T.push_back(temp);
	}

	//Pb_test and Tb_test
	if (basis_type_test==ONE_DIM_LINEAR)//baisi function type
	{		
		Pb_test = P;
		Tb_test = T;
		number_of_local_basis_test = 2;	
	}
	else if (basis_type_test == ONE_DIM_QUADRATIC)
	{
		Pb_test.push_back(vector<double>{x_min});
		for (int i = 0; i < 2*number_of_elements; ++i)
		{
			Pb_test.push_back(vector<double>{(x_max - x_min) / (2*number_of_elements) * (i + 1)});
			vector<int> temp;		
		}
		for (int i=0;i<number_of_elements;++i)
		{
			vector<int> temp;
			temp.push_back(2 * (i + 1) - 1 - 1);
			temp.push_back(2 * (i + 1) + 1 - 1);
			temp.push_back(2 * (i + 1) - 1);
			Tb_test.push_back(temp);
		}

		number_of_local_basis_test = 3;
	}
	
	else
	{

	}
	//Pb_trial and Tb_trial
	if (basis_type_trial == ONE_DIM_LINEAR)//baisi function type
	{
		
		Pb_trial = P;
		Tb_trial = T;
		number_of_local_basis_trial = 2;
	}
	else if (basis_type_trial == ONE_DIM_QUADRATIC)
	{
		Pb_trial.push_back(vector<double>{x_min});
		for (int i = 0; i < 2 * number_of_elements; ++i)
		{
			Pb_trial.push_back(vector<double>{(x_max - x_min) / (2 * number_of_elements) * (i + 1)});
			vector<int> temp;
		}
		for (int i = 0; i < number_of_elements; ++i)
		{
			vector<int> temp;
			temp.push_back(2 * (i + 1) - 1 - 1);
			temp.push_back(2 * (i + 1) + 1 - 1);
			temp.push_back(2 * (i + 1) - 1);
			Tb_trial.push_back(temp);
		}

		number_of_local_basis_trial = 3;
	}
	else
	{

	}
	
	return true;	
}

//generate the boundary conditions contains the 
vector<boundary_condition> generate_boundarynodes(const int ele_num)
{
	vector<boundary_condition> boundarynodes;
	boundary_condition temp;

	//example 1 
	//boundary for the left side
	/*temp.type = Dirichlet;
	temp.glo_node_index = 0;
	temp.value = vector <double>{ 0 };
	boundarynodes.push_back(temp);

	
	//boundary for the right side of
	temp.type = Dirichlet;
	temp.glo_node_index = ele_num;
	temp.value = vector <double>{ cos(1.0) };
	boundarynodes.push_back(temp);*/

	//example 2 
	//boundary for the left side
	temp.type = Dirichlet;
	temp.glo_node_index = 0;
	temp.value = vector <double>{ 0 };
	boundarynodes.push_back(temp);


	//boundary for the right side of
	temp.type = Dirichlet;
	temp.glo_node_index = 2*ele_num;
	temp.value = vector <double>{ cos(1.0) };
	boundarynodes.push_back(temp);

	//example 3
	//boundary for the left side
	/*temp.type = Dirichlet;
	temp.glo_node_index = 0;
	temp.value = vector <double>{ 0 };
	boundarynodes.push_back(temp);

	//boundary for the right side
	temp.type = Neumann;
	temp.glo_node_index = 2 * ele_num;
	temp.value = vector<double>{ cos(1.0) - sin(1) };
	temp.outer_normal_vector = vector<double>{ 1.0,0.0,0.0 };//outer normal vector of the right side boundary
	boundarynodes.push_back(temp);*/


	//example 4
	//boundary for the left side
	/*temp.type = Robin;
	temp.glo_node_index = 0;
	temp.value = vector <double>{1,1};
	temp.outer_normal_vector = vector<double>{ -1.0,0.0,0.0 };//outer normal vector of the left side boundary
	boundarynodes.push_back(temp);

	//boundary for the right side
	temp.type = Dirichlet;
	temp.glo_node_index = 2 * ele_num;
	temp.value = vector<double>{ cos(1.0)};
	boundarynodes.push_back(temp);*/

	return boundarynodes;
}

//local basis function for one dimension
//basis_type: define the basis type
//vertices: the coordinates of the vertices
//basis_index: index of the basis function, first one and second one for 1D problem
//basis_der_x: derivative order of the basis function for x
//The definition of basis local function for both test and trial

double FE_basis_local_func_1D(vector<double> x,vector<vector<double>> vertices, BASIS_TYPE basis_type,int basis_index, int basis_der_x)
{
	double h = vertices[1][0] - vertices[0][0];
	if (basis_type ==ONE_DIM_LINEAR)//different basis type 
	{
		if (basis_index == 0)//the first basis function
		{
			if (basis_der_x==0)//zero derivative of the basis function
			{
				return (vertices[1][0] - x[0]) / h;
			}
			else if(basis_der_x ==1)//first derivation of the basis function
			{
				return -1 / h;
			}
			else if (basis_der_x >= 2)
			{
				return 0;
			}
			else
			{
				cout << "\n error: wrong input for basis derivative\n";
			}
		}
		else if (basis_index == 1)//the second basis function
		{
			if (basis_der_x == 0)
			{
				return (x[0]- vertices[0][0]) / h;
			}
			else if (basis_der_x == 1)
			{
				return 1 / h;
			}
			else if (basis_der_x >= 2)
			{
				return 0;
			}
			else
			{
				cout << "\n error: wrong input for basis derivative\n";
			}
		}
		else
			cout << "error: wrong input for basis index\n";

	}
	else if (basis_type == ONE_DIM_QUADRATIC)
	{
		if (basis_index == 0)//the first basis function
		{
			if (basis_der_x == 0)//zero derivative of the basis function
			{
				return 2 * pow((x[0] - vertices[0][0]) / h, 2.0) - 3 * (x[0] - vertices[0][0]) / h + 1;
					
			}
			else if (basis_der_x == 1)//first derivation of the basis function
			{
				return 4*(x[0]- vertices[0][0])/h/h-3.0/h;
			}
			else if (basis_der_x == 2)
			{
				return 4/h/h;
			}
			else if (basis_der_x >=3)
			{
				return 0;
			}
			else
			{
				cout << "\n error: wrong input for basis derivative\n";
			}
		}
		else if (basis_index == 1)//the second basis function
		{
			if (basis_der_x == 0)
			{
				return 2*pow((x[0] - vertices[0][0]) / h,2)- (x[0] - vertices[0][0])/h;
			}
			else if (basis_der_x == 1)
			{
				return 4 * (x[0] - vertices[0][0]) / h / h-1/h;
			}
			else if (basis_der_x == 2)
			{
				return 4/h/h;
			}
			else
			{
				cout << "\n error: wrong input for basis derivative\n";
			}
		}
		else if (basis_index==2)//the third basis function
		{
			if (basis_der_x == 0)
			{
				return -4*pow((x[0] - vertices[0][0]) / h,2)+4*(x[0] - vertices[0][0])/h;
			}
			else if (basis_der_x == 1)
			{
				return -8* (x[0] - vertices[0][0]) / h / h +4/h;
			}
			else if (basis_der_x == 2)
			{
				return -8/h/h;
			}
			else if (basis_der_x >= 3)
			{
				return 0;
			}

			else
			{
				cout << "\n error: wrong input for basis derivative\n";
			}
		}
		else
			cout << "error: wrong input for basis index\n";

	}
	return 0;
}

double Gauss_quadrature_1D_trial_test(double (*coe)(vector<double>),vector<double> Gauss_weights,vector<vector<double>> Gauss_nodes,vector<vector<double>> vertices, 
	BASIS_TYPE basis_type_trial, int basis_index_trial, int basis_der_x_trial,
	BASIS_TYPE basis_type_test, int basis_index_test, int basis_der_x_test)
{
	int Gpn = Gauss_weights.size();
	double int_value = 0;
	for (int k = 0; k < Gpn; ++k)
	{
		int_value += Gauss_weights[k]* coe(Gauss_nodes[k])*FE_basis_local_func_1D(Gauss_nodes[k], vertices, basis_type_trial, basis_index_trial, basis_der_x_trial)*FE_basis_local_func_1D(Gauss_nodes[k], vertices, basis_type_test, basis_index_test, basis_der_x_test);
	}
	return int_value;
}
double Gauss_quadrature_1D_load_trial(double(*load)(vector<double>), vector<double> Gauss_weights, vector<vector<double>> Gauss_nodes, vector<vector<double>> vertices,
	BASIS_TYPE basis_type_test, int basis_index_test, int basis_der_x_test)
{
	int Gpn = Gauss_weights.size();
	double int_value = 0;
	for (int k = 0; k < Gpn; ++k)
	{
		int_value += Gauss_weights[k] * load(Gauss_nodes[k])*FE_basis_local_func_1D(Gauss_nodes[k], vertices, basis_type_test, basis_index_test, basis_der_x_test);
	}
	return int_value;
}

//coefficients
double coe_function(vector<double> x)
{
	return exp(x[0]);
}
double load_function(vector<double> x)
{
	return -exp(x[0])*(cos(x[0]) - 2 * sin(x[0]) - x[0] * cos(x[0]) - x[0] * sin(x[0]));
}

void generate_Gauss_local(vector<vector<double>> vertices, int Gauss_type, vector<vector<double>> &Gauss_nodes, vector<double> &Gauss_weights)
{
	//Gauss quadrature points on [-1,1]
	vector<double> w;
	vector<double> p;
	if (Gauss_type == 2)
	{
		w = vector<double>{1,1};
		p = vector<double>{ -1 / sqrt(3.0), 1 / sqrt(3.0) };
	}
	
	else if (Gauss_type == 4)
	{
		w = vector<double>{ 0.3478548451, 0.3478548451, 0.6521451549, 0.6521451549 };
		p = vector<double>{ 0.8611363116, -0.8611363116, 0.3399810436, -0.3399810436 };
	}
		
	else if (Gauss_type == 8)
	{
		w = vector<double>{ 0.1012285363, 0.1012285363, 0.2223810345, 0.2223810345, 0.3137066459, 0.3137066459, 0.3626837834, 0.3626837834 };
		p = vector<double>{ 0.9602898565, -0.9602898565, 0.7966664774, -0.7966664774, 0.5255324099, -0.5255324099, 0.1834346425, -0.1834346425 };
	}
	else
	{
		cout << "\n error: Gauss quadrature type error\n";
		return;
	}

	//Gauss quadrature point on [vertices[0],vertices[1]]
	for (unsigned int i=0;i<w.size();++i)
	{
		if (vertices[0].size()==1)//one dimension
		{
			Gauss_weights.push_back((vertices[1][0] - vertices[0][0]) / 2.0*w[i]);
			Gauss_nodes.push_back(vector<double> {(vertices[1][0] - vertices[0][0]) / 2.0*p[i] + (vertices[1][0] + vertices[0][0]) / 2.0});

		}
		else if (vertices[0].size() == 2)//two dimension
		{

		}
		else if (vertices[0].size() == 3)//three dimension
		{

		}
		else
		{
			cout << "\nerror: element node coordinates error\n";
		}
	}
	
}

SpMat assemble_matrix_1D(pair<int,int>matrix_size,int number_of_elements,vector<vector<double>> P, vector<vector<int>>T, vector<vector<int>>Tb_trial, vector<vector<int>>Tb_test,int number_of_local_func_trial, int number_of_local_func_test,
	BASIS_TYPE basis_type_trial, int basis_der_x_trial, BASIS_TYPE basis_type_test, int basis_der_x_test,int Gauss_type)
{
	SpMat A(matrix_size.first, matrix_size.second);
	
	for (int n = 0; n < number_of_elements; ++n)
	{
		vector<vector<double>> vertices;//vertices (coordinates) of nth element

		for (unsigned int i=0;i<T[n].size();++i)
		{
			vertices.push_back(P[T[n][i]]);
		}
		vector<double> Gauss_weights;
		vector<vector<double>> Gauss_nodes;
		//generate the Gauss integration points
		generate_Gauss_local(vertices, Gauss_type,  Gauss_nodes, Gauss_weights);
		for (int alpha=0;alpha<number_of_local_func_trial;++alpha)
			for (int beta=0;beta<number_of_local_func_test;++beta)
			{
				double int_value = Gauss_quadrature_1D_trial_test(coe_function,Gauss_weights, Gauss_nodes, vertices,
					basis_type_trial, alpha,  basis_der_x_trial,
					 basis_type_test,  beta, basis_der_x_test);
				int i = Tb_test[n][beta];
				int j = Tb_trial[n][alpha];
				A.coeffRef(i, j) += int_value;
			}
	}
	return A;
}

Rsv assemble_vector_1D(pair<int, int>matrix_size,int number_of_elements, vector<vector<double>> P, vector<vector<int>>T, vector<vector<int>>Tb_test, int number_of_local_func_test, BASIS_TYPE basis_type_test, int basis_der_x_test, int Gauss_type)
{
	auto b = Rsv(matrix_size.first);

	for (int i=0;i< matrix_size.first;++i)
	{
		b(i) = 0;
	}

	for (int n = 0; n < number_of_elements; ++n)
	{
		vector<vector<double>> vertices;//vertices of nth element

		for (unsigned int i = 0; i < T[n].size(); ++i)
		{
			vertices.push_back(P[T[n][i]]);
		}
		vector<double> Gauss_weights;
		vector<vector<double>> Gauss_nodes;
		//generate the Gauss integration points
		generate_Gauss_local(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		for (int beta = 0; beta < number_of_local_func_test; ++beta)
		{
			double int_value = Gauss_quadrature_1D_load_trial(load_function, Gauss_weights, Gauss_nodes, vertices, basis_type_test, beta, basis_der_x_test);

			b[Tb_test[n][beta]] += int_value;
		}
	}
	return b;
}

void treat_boundary(SpMat &A, Rsv &b, vector<boundary_condition> boundarynodes,vector<vector<double>> Pb_test)
{
	for (unsigned int k=0;k<boundarynodes.size();++k)
	{
		if (boundarynodes[k].type==Dirichlet)
		{
			int i = boundarynodes[k].glo_node_index;
			for (int j=0;j<A.cols();++j)
			{
				A.coeffRef(i, j) = 0;
			}
			A.coeffRef(i, i) = 1;

			b[i] = boundarynodes[k].value[0];
		}
		else if (boundarynodes[k].type == Neumann)
		{
			int i = boundarynodes[k].glo_node_index;
			b[i] += boundarynodes[k].value[0]*coe_function(Pb_test[i])*boundarynodes[k].outer_normal_vector[0];
		}
		else if (boundarynodes[k].type == Robin)//corresponding to the outer normal vector of the boundary
		{
			int i = boundarynodes[k].glo_node_index;
			b[i] += boundarynodes[k].value[1] * coe_function(Pb_test[i])*boundarynodes[k].outer_normal_vector[0];

			A.coeffRef(i, i) += boundarynodes[k].value[0] * coe_function(Pb_test[i])*boundarynodes[k].outer_normal_vector[0];
		}
		else
		{

		}
	}

}


double local_FE_function_1D(vector<double>uh_vector,int number_of_local_func, vector<double>x,vector<vector<double>> vertices, BASIS_TYPE basis_type, int basis_der_x)
{
	double result = 0;
	for (int k = 0; k < number_of_local_func; ++k)
	{
		result += uh_vector[k] * FE_basis_local_func_1D(x, vertices, basis_type, k, basis_der_x);

	}
	return result;
}



//Compute the error with Gauss integration
double Gauss_int_error_1D(double(*ana_sol)(double,int),vector<double> uh_vector,vector<double> Gauss_weights, vector<vector<double>> Gauss_nodes, vector<vector<double>> vertices,
	BASIS_TYPE basis_type_trial, int s, int number_of_local_func)
{
	int Gpn = Gauss_weights.size();
	double int_value = 0;
	for (int k = 0; k < Gpn; ++k)
	{
		int_value += Gauss_weights[k] * pow(ana_sol(Gauss_nodes[k][0],s)- local_FE_function_1D(uh_vector, number_of_local_func, Gauss_nodes[k], vertices, basis_type_trial,s),2);
	}
	return int_value;
}

double analytical_solution(double x, int s)
{
	if (s == 0)
	{
		return x * cos(x);

	}
	else if (s == 1)
	{
		return cos(x) - x * sin(x);
	}
}

//Compute the L2 and H1 error
//s=0: L2
//s=1: H1
double compute_Hs_error(int number_of_elements, vector<vector<double>> P, vector<vector<int>>T, vector<vector<int>>Tb_test, 
	int number_of_local_func_trial, BASIS_TYPE basis_type_trial,int s, int Gauss_type,vector<double> solution_vector)
{
	double error = 0;
	for (int n = 0; n < number_of_elements; ++n)
	{
		vector<vector<double>> vertices;//vertices (coordinates) of nth element

		for (unsigned int i = 0; i < T[n].size(); ++i)
		{
			vertices.push_back(P[T[n][i]]);
		}
		vector<double> Gauss_weights;
		vector<vector<double>> Gauss_nodes;
		//generate the Gauss integration points
		generate_Gauss_local(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		
		//local solution vector
		vector<double> uh_vector;
		for (int i=0;i< Tb_test[n].size();++i)
		{
			uh_vector.push_back(solution_vector[Tb_test[n][i]]);
		}
		
		double int_value = Gauss_int_error_1D(analytical_solution, uh_vector,Gauss_weights, Gauss_nodes, vertices, basis_type_trial, s, number_of_local_func_trial);;

		error += int_value;
	}
	return sqrt(error);
}

double compute_L_infi_error(int number_of_elements, vector<vector<double>> P, vector<vector<int>>T, vector<vector<int>>Tb_test,
	int number_of_local_func_trial, BASIS_TYPE basis_type_trial, int basis_der_x,int Gauss_type, vector<double> solution_vector)
{
	
	vector<double> r_vec;
	for (int n = 0; n < number_of_elements; ++n)
	{
		vector<vector<double>> vertices;//vertices (coordinates) of nth element

		for (unsigned int i = 0; i < T[n].size(); ++i)
		{
			vertices.push_back(P[T[n][i]]);
		}

		vector<double> Gauss_weights;
		vector<vector<double>> Gauss_nodes;
		//generate the Gauss integration points
		generate_Gauss_local(vertices, Gauss_type, Gauss_nodes, Gauss_weights);

		//local solution vector
		vector<double> uh_vector;
		for (int i = 0; i < Tb_test[n].size(); ++i)
		{
			uh_vector.push_back(solution_vector[Tb_test[n][i]]);
		}

		//initial value is the first Gauss_nodes
		
		double r = analytical_solution(Gauss_nodes[0][0], 0) - local_FE_function_1D(uh_vector, number_of_local_func_trial, Gauss_nodes[0], vertices, basis_type_trial, basis_der_x);
		for (int k=0;k< Gauss_nodes.size();++k)
		{
			double dr = fabs(analytical_solution(Gauss_nodes[k][0], 0) - local_FE_function_1D(uh_vector, number_of_local_func_trial, Gauss_nodes[k], vertices, basis_type_trial, basis_der_x));
			if (r < dr)
			{
				r = dr;
			}
		}
		r_vec.push_back(r);

	}
	return *max_element(r_vec.begin(), r_vec.end());
}




int main()
{
	int ele_num = 4;
	const double x_min = 0;
	const double x_max = 1;
	const BASIS_TYPE basis_type_trial= ONE_DIM_QUADRATIC;
	const BASIS_TYPE basis_type_test = ONE_DIM_QUADRATIC;

	const int basis_der_x_trial = 1;
	const int basis_der_x_test = 1;

	//Gauss quadrature type is 4
	const int Gauss_type = 4;

	cout << "\nFirst homework of the finite element method course: 1D elliptic equation with linear basis function \n";
	cout << "\nSimulation begins\n";
	cout << "\nInput the number of the elements:\n";
	cin >> ele_num;

	vector<vector<double>> P;
	vector<vector<int>> T;
	vector<vector<double>> Pb_trial; vector<vector<int>> Tb_trial; int number_of_local_basis_trial;
	vector<vector<double>> Pb_test; vector<vector<int>> Tb_test; int number_of_local_basis_test;

	//Generate the P, T matrix for mesh, finite element of trial and test
	 mesh_generate(basis_type_trial,basis_type_test,x_min,x_max, ele_num,P,T,Pb_trial,Tb_trial,number_of_local_basis_trial,Pb_test,Tb_test,number_of_local_basis_test);

	//Generate the boundarynodes
	vector<boundary_condition> boundarynodes = generate_boundarynodes(ele_num);

	pair<int, int> matrix_size{Pb_test.size(),Pb_trial.size()};	
	
	//Assemble matrix
	SpMat A= assemble_matrix_1D(matrix_size, ele_num, P, T, Tb_trial, Tb_test, number_of_local_basis_trial, number_of_local_basis_test, basis_type_trial, basis_der_x_trial, basis_type_test, basis_der_x_test, Gauss_type);
	
#ifdef _DEBUG
	for (int i=0;i<A.rows();++i)
	{
		for (int j=0;j<A.cols();++j)
		{
			cout << setw(10)<<A.coeffRef(i, j)<<"\t";
		}
		cout << endl;

	}
#endif

	//Assemble right side vector
	Rsv b = assemble_vector_1D(matrix_size, ele_num, P, T, Tb_test, number_of_local_basis_test, basis_type_test, 0, Gauss_type);
#ifdef _DEBUG
	for (int i=0;i<b.size();++i)
	{
		cout << b[i] << endl;
	}

#endif
	//treat boundary conditions
	treat_boundary(A, b, boundarynodes,Pb_test);

	//Solve the linear equations with Eigen
	//Define a SparseLU solver
	Eigen::SparseLU<SpMat> solver;
	solver.compute(A);
	if (solver.info() != Eigen::Success)
	{
		cout << "\nerror: decomposition failed when using Eigen to compute matrix A\n";
	}
	Rsv x = solver.solve(b);
	if (solver.info() != Eigen::Success)
	{
		cout << "\nerror: Solving failed when using Eigen to solve Ax=b\n";
	}	

	//Print the solution
	cout << "\nSolution:\n";
	for (int i=0;i<x.size();++i)
	{
		cout << "x[" << i << "]=\t" << setiosflags(ios::scientific)<<x(i) << endl;
	}

	//Print the error of the solution
	cout << "\nError of all nodes:\n";
	double max_error = fabs(x(0) - analytical_solution(P[0][0],0));
	for (int i = 0; i < x.size(); ++i)
	{
		
		cout << "x["<<i<<"]=" << setiosflags(ios::scientific)<< fabs(x(i)- analytical_solution(Pb_trial[i][0],0))<< endl;
		if (max_error < fabs(x(i) - analytical_solution(Pb_trial[i][0],0)))
		{
			max_error = fabs(x(i) - analytical_solution(Pb_trial[i][0],0));
		}
	}
	cout << "The maximum absolute error at all nodes:\t" << setiosflags(ios::scientific) << max_error << endl;


	vector <double> solution_vec;
	for (int i = 0; i < x.size(); ++i)
	{
		solution_vec.push_back(x[i]);
	}

	//Compute the L inifite norm erro
	cout << "L∞ error=\t" << compute_L_infi_error(ele_num, P, T, Tb_trial, number_of_local_basis_trial, basis_type_trial, 0, Gauss_type, solution_vec) << endl;

	//Compute the L2 norm error and H1 semi-norm error
	
	cout << "L2 error=\t"<<compute_Hs_error(ele_num, P, T, Tb_trial, number_of_local_basis_trial, basis_type_trial, 0, Gauss_type, solution_vec)<<endl;
	cout << "H1 error=\t" << compute_Hs_error(ele_num, P, T, Tb_trial, number_of_local_basis_trial, basis_type_trial, 1, Gauss_type, solution_vec) << endl;

	//Output the solution to files and plot the solution with plot.py	
	string filename = string("solution_") + string("1_") + to_string(ele_num) + string(".txt");
	ofstream f(filename);
	f << "h=" << (x_max - x_min) / ele_num<<"\t max_error="<<max_error<<endl;
	f << "x\t" << "analytical solution\t"<<"finite element solution" << endl;
	for (int i=0;i<x.size();++i)
	{
		f << Pb_trial[i][0] << "\t" << analytical_solution(Pb_trial[i][0],0) << "\t" << x(i) << endl;
	}
	f.close();
	return 0;
}
